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Abstract 



We study the dynamics of a stepped crystal surface during evaporation, 

using the classical model of Burton, Cabrera and Frank, in which the dynamics 

of the surface is represented as a motion of parallel, monoatomic steps. The 

validity of the continuum approximation treated by Frank is checked against 

numerical calculations and simple, qualitative arguments. The continuum 

approximation is found to suffer from limitations related, in particular, to the 

existence of angular points. These limitations are often related to an adatom 

detachment rate of adatoms which is higher on the lower side of each step 

than on the upper side ("Schwoebel effect"). 
Pacs numbers: 68.45.Da; 05.70.Ln; 68.20 
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I. INTRODUCTION 



The study of the dynamics of a crystal surface under nonequilibrium conditions is an 
old subject, the interest of which was renewed in the last years due to modern experimen- 
tal techniques, such as electron reflection microscopy [|TJ @, which allow a more detailed 
observation of dynamic phenomena of the surfaces both during evaporation and growth. 
Moreover, equilibrium is practically never reached, and even a surface which is apparently 
at equilibrium on a certain lengthscale, will present nonequilibrium shapes at larger length- 
scales. Therefore, it is of interest to understand how these features can be formed, either 
during growth or during annealing. 

Burton, Cabrera and Frank || (BCF) have developed a simple theory to describe the 
growth of a stepped, dislocation-free crystal surface. They wrote equations of motion for 
the steps. The purpose of the present article is to solve these equations for special cases 
corresponding to a given initial profile of the surface. 

When trying to solve the BCF equations, a possible approach is to make the continuum 
approximation. Then, as will be seen in Section 2, it is possible to use a theorem discovered 
by Frank [|J. However, the continuum approximation is questionable if the slope of the 
surface is discontinuous. For that reason, the present work is focussed on the evaporation 
of surfaces which initially contain corners (Fig. [I] a). 

Only evaporation will be considered here, in order to avoid problems due to nucleation 
of new terraces during growth. One can also notice that in situ study of surface dynamics 
is more easily performed under evaporation than under growth conditions Q @]. In this 
paper we investigate, within the BCF theory, the evaporation of a defect-free surface made 
of parallel, but not equidistant steps. Our main motivation is to test, in particular instances, 
the validity of the continuum approximation developed by Frank Q. The continuum ap- 
proximation will be shown to fail when curvature changes abruptly at some place on the 
surface. In the original BCF paper, atoms detaching from a step were assumed to be emitted 
with equal probability from the upper and from the lower terrace. This will be called the 
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symmetric case. However, as noticed by Ehrlich ||, a step may be partially rather than 
totally absorbing. Generally, one side is more absorbing than the other |J. This is known 
as Schwoebel effect. 

We will consider solid surfaces made of parallel steps. The structure of the surface at time 
t is fully characterized by the positions x n (t) of the successive steps labelled by the index 
n = 1, 2, ... as shown in Fig. [I] a. The functions x n (t) satisfy the following equation, which 
follows from the BCF theory, appropriately modified to take into account the Schwoebel 
effect: 

X n = -(fi(x n - X n -i) - <p r (x n+ i - X n ) (1) 

The expressions of ipi(l) and <p r (l) are derived in Reference and, for completeness, 

briefly rederived in Appendix A. The first one, 

cosh.(nl) - 1 + jysmh(Kl) 
(l + sinh(KZ) + + f£) cosh(/d) ' 

gives the net flux of outgoing adatoms from a step to the upper terrace, and 

cosh(ft/) — 1 + j^sinh(«Z) 
(l + §gn) sinh(KZ) + (*p + !§) cosh(/d) ' 

is the net flux of outgoing adatoms from a step to the lower terrace. In these expressions, 
Po is the equilibrium density of adatoms on the high symmetry surface at the appropriate 
temperature and vapour pressure, D is the surface diffusion constant of adatoms and 1/Vq 
is the evaporation probability of an adatom per unit time. F is zero in the cases addressed 
here, but it would be the beam intensity in the case of growth by molecular beam epitaxy. 
D' is the unit time probability of an adatom to stick to a step when it is just aside the step, 
and D" is the unit time probability of an adatom to stick to a step when it is just above that 
step (Fig. |T] b). The situation where D' ^ D and D" ^ D will correspond to situations in 
which the sticking probability of adatoms to steps will depend on the terrace towards which 
they move. More specifically, we will call normal Schwoebel effect when particles mainly 
move towards the upper terrace, D" << D, while the situation in which particles move 
towards the lower terrace will be named inverse Schwoebel effect. Finally, 



<Pi{l) = (Po - Ft ) j- ■ — . — — kD (2) 



Ml) = (Po - Ft ) , ■ „ n2X ■ f- ■ —r — — kD (3) 



is the reciprocal of the average distance on which an adatom would diffuse on an infinite 
terrace before evaporating. This distance is always much larger than the interatomic spacing, 
taken to be the length unit throughout this work, while the time an adatom needs to diffuse 
along the interatomic distance is the unit time. Normal Schwoebel effect is known to produce 
instabilities during evaporation , but the present study is restricted to cases where no such 
instability occurs. 

As seen from (0) and (|^), tpi(l) and v 2 r(0 are proportional to / for small I and constant 
for large I. In the symmetric case one has D' = D" = D and, since k << 1, (0) and ([3]) 
reduce to the BCF formula 



The conditions of applicability of the BCF theory, which implies the stability of the step 
flow regime, are the following. 

1. No dislocations should be present. Generally this implies that only a small part of the 
surface is observed, e.g. 0.01 cm x 0.01 cm for Si wafers 0. 

2. No surface vacancies should be formed on terraces. This implies that the temperature 
should not exceed a certain threshold, e.g. 1200 K for Si(001) M though much higher 



3. The vicinal orientation should be stable with respect to facet formation. 

4. The steps should not undergo instabilities such as the one discovered by Bales and 
Zangwill f7| in the case of growth. 

In the next section, the continuum approximation of ([[]) will be introduced and exploited. 
We will then focus our attention on the evaporation of a surface limited initially by two semi- 
infinite half planes as sketched in Fig. [I] a. One of them is assumed to be a high symmetry 





for Si(lll) [0]. 



surface and the other one is made of parallel steps. In section 3, the long time behaviour of 
the solution of ([I]) will be derived analytically. In section 4, exact solutions will be derived 
for approximations of ([!]). Numerical solutions of (|1|) will be presented in section 5 and will 
be compared with the continuum approximation. Finally, more complicated initial profiles 
will be investigated in Section 6, namely a periodic surface with grooves. In the discussion 
we summarize our main results, and technical details are worked out in the appendices. 



II. THE CONTINUUM MODEL AND FRANK CONSTRUCTION 

In the continuum approximation, the step position x n is regarded as a continuous and 
generally different iable function of the local surface height — n, also considered as a continu- 
ous variable z. Therefore, x n — > x(—z). Moreover, introducing the derivatives x' = —dx/dz, 
x" = d 2 x/dz 2 , etc., and <p'(x') = dip(x')/dx', and assuming x'" and the further derivatives 
to be small, one can write (|1|) as 

T " T iii r "2 

x = -24>{x') + X - - —${x') - X —<P"{x>) (6) 

where 

m = jgK z xQ+jg : (-sO (7) 

i;{x') = ipti-x') - <pr{-a/) (8) 

and tp' and ip" are the derivatives of ip. 

In this section, we consider the approximation obtained by keeping only the first term 
in®, 

x = -2(j){x') (9) 

This approximation will be called the continuum approximation. It clearly fails near a corner, 
when x" is infinite. Moreover, it may be expected to be better when there is no Schwoebel 
effect, since then the second term of the right hand side of eq.(O) vanishes according to (R). 
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Equation @ may alternatively be written as a relation between z' = (5z/5x)t = 
l/(dx/dz) t = —x' and z = (dx/dz) x = x/x': 

z = 2*>(i) (10) 
z' 

From this expression one can deduce the decay rate of the surface, v, at a point charac- 



terized by a direction vector n = (n x ,0,n z ), with n x = —1/yl + x' 2 . Due to the geometry, 
the decay rate is given by v = z n z = x n x , so that 

= _ 2 = /(n) (11) 

This equation shows that in the continuum limit the decay rate of the surface depends 
only on its local orientation. The problem of the decay or growth of a surface when its 
velocity is a function of the orientation has been investigated by Frank ||. Frank proved 
that the points at the surface with a given normal orientation move on straight lines. Let 
us consider the initial profile of Fig. [I] a. For a velocity given by eq. (|TT|) ? we can draw the 



polar plot T of the inverse of the velocity, as is shown in Fig. |2|. From the curve (r) we can 
deduce the further evolution of the profile. If we take a point P of the crystal where the 
normal is n, it determines a point M of (r) through OM / /ft. Then Frank's theorem states 
that the point P moves on a straight line parallel to the normal n' to (r) at M. We will 
come back to this property later on. 

In order to test the validity of the continuum approximation, let us consider the initial 
profile of Fig. || a such that all initial terraces have the same width except the first one. If 
there is a total normal Schwoebel effect such that there is no interaction between a step and 
its upper terrace, that is ipi(l) = 0, then all steps will move at the same constant velocity, 
and therefore terrace widths will remain constant. In other case, the interaction of the 
second step with a different one will induce a modification in the step velocity, and then 
the velocity of the steps will change in time and will be different from one to another. The 
continuum approximation is not sensitive to this qualitatively different behaviour induced 
by the Schwoebel effect. From equation (0), one sees that the function <p(x) does not change 
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qualitatively whether <pi(—x') is zero or not. This effect will be studied quantitatively in the 
next section. 



III. UPPER LEDGES AT LONG TIMES 

We are now concerned with the evolution of the simple profile already considered in 
section 2. At the beginning, the distances li,l2,--- are finite and equal. It will be shown 
self-consistently that the distances go to infinity for long times. We take this as an Ansatz 
to be proved later. For large values of I, equations (0) and (|^) read 

(pi(l) = A(A' - 2e- Kl + C'e- 2Kl ) (12) 
<p r (l) = A (A" - 2e~ Kl + C"e~ 2Kl ) (13) 

where A, A', A" , C and C" are given in appendix A. It is necessary to go to second order 
in exp(— kI) because, when D" goes to zero, A" and C" go to infinity while A goes to zero. 
A similar effect takes place when D' = 0. 
From (P, one obtains 

l n = X n+ i - X n = -(pi(l n ) - (fr(ln+l) + <Pl(L-l) + <Pr(L) (14) 

Combining eq.(Tl) with eqs. ( P^|) and (0), one obtains for n > 1 



l n = A (2e" K/ " +1 - 2e" K/ "- 1 - C'e~ 2Kln + C'e' 2 ^- 1 - C" e - 2Kln+1 + C"'e" 2K '") (15) 
and for n = 1: 

h = A (2e~ Kh + (C - C')e- 2Kh - C"e' 2Kh ) (16) 

If C and C" are finite, the terms in exp(— 2nl) may be neglected in (|16D for large I. Then 
these equations turn out to have solutions of the form 

l n (t) = -ln(B n t) (17) 

K 



Indeed, insertion of (|17D into Qlq) yields 



B n+ i B n _i 2An 
while relations ([T7|) and flIE|) yield 

1 1 



B 2 2Ak < 19 » 



From ( |18|) and fll9D one deduces 

1 n 



(20) 



5 2n 2kA 

1 1 n 

ir + 7T-r (21) 



5 2n+1 fix 2kA 

The numerical solution (section 5) shows that the correct solution is the most symmetrical 
one, namely Bi = 4kA, so that 

l n (t) -lln(-) (22) 



K \n 



Under total normal Schwoebel effect, D" = 0, then ipi = 0, as seen from (|2|). Then it 
follows from fll4] ) that l n = —ipi(l n+ i) + <fi(l n ) and, if all Z n 's are equal at the beginning, 
they remain equal, as expected. If total inverse Schwoebel effect is considered, D' = 0, then 
cp r = as seen from (||), and for kI » 1 and kD < D": 

(ptil) = (p - Ft ) tanhfcZ) ~ (p - Fr )(l - 2e Kl ) (23) 

In this case, expression ( |T4"D reduces to 

L = -Pliln) + Vliln-l) (24) 

which has a solution similar to eq.(|22|) 

l n (t) = ±ln(-) (25) 



2k \n 

These logarithmic solutions correspond to a self-similar shape of the surface, since l an (at) 
l n (t) for any value of the parameter a. 
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IV. PIECEWISE LINEAR APPROXIMATION 



A. The approximation 

In this section we are interested in the evolution of the same profile as the one considered 
in the previous section, but we will focus our attention on the short time evolution. We are 
again interested in solving eqs.(P subject to the initial conditions (Fig. |l] a) 

Xn - x n -i = K Q ) ,n> 1 

x\ — xo = oo (26) 

with 1(0) < 2/k, so that the saturation regime has not been reached. Otherwise, the result 
of the previous section would apply. 

In order to solve the evolution equations, we will take advantage of the form of the 
hyperbolic tangent, which is approximated by the piecewise linear function 



We will study the evolution of the initial profile eq.pq) within the piecewise linear 
approximation eq.(p7|) in both the symmetric model and under total inverse Schwoebel 
effect. 



B. Symmetric model 

In this case, when all Z n 's are smaller than 1/k, using the form of our initial conditions 
(0), eqs.(|) read 

x x = -—^-(x 2 - xi) - Ak 

i n = -^-(x n+1 - X n -i) ,n>l (28) 

where A = D(po — Ftq). As usual, it is convenient to rewrite the evolution equations in 
terms of the step widths, 
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Tl > 1 



(29) 



The solution of eqs. (p9|) , as derived in appendix B, is of the form 



kn(t) 



n-l 



Jo(AkH) + J 2n (AK 2 t) + MAkH) 



1=1 



K 



1(0) 



hn+l{t) = h(AK 2 t) - J X (AkH) + J 2n +l(AK 2 t) + 



n-l 



+ J2 J 2i+l(AK 2 t) 



i=l 



--1(0) 



h(t) = l(0) + 



K 

oo 



J 1 (AK 2 t)+2j2J2 i+ i(AK 2 t) 
1=1 



n > 1 



-1(0) 



(30) 



This solution shows that at the beginning the first step moves faster than the others, 
and the rest of the steps move faster progressively. 

However, the validity of eq.(|30"D is restricted by the fact that in eq.(p8|) we have assumed 
all / ra 's to be smaller than \ j k . Thus, when the width of one step surpasses the characteristic 
length 2/k, this solution will break down. The first step to fulfill this requirement is the first 
step, h, for a time t ~ 2.Q5/Ak 2 . At longer times, ( pO] ) does not apply. 



C. Total inverse Schwoebel effect 



The piecewise linear approximation will now be applied to the case (p r — within the 



piecewise linear approximation. Taking the initial conditions given by eq.([26|), as kD « D", 
the evolution equations read, at least for small times, 



x\ = —Ak 

x n Ak (x n x x ~i) 



n>2 



(31) 



In terms of the widths, one can rewrite these equations as 



i x = Ak(1 - nli) 

i n = Ak 2 ^ - i r 



,n>2 



(32) 
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which, using as initial conditions eqs.fl26|), has a solution of the form 



Inif) = - 
K 



K 



-An 2 t 



e £^ (33) 

p=l t>- 



In contrast with eq. (|30|) , this solution turns out to be valid for all times, because l n {t) < 
1/ ' k for all t and n, and therefore the regime of the previous section, valid in the symmetric 
model, is never achieved in the completely asymmetric one. However, its behaviour is 
qualitatively similar to the symmetric model. At short times the first step begins to move 
faster, while the others move at the same speed. After some time, the second step increases 
its velocity, separates from the following steps and approximates the first step, and so on. 
This behaviour is easily seen studying the difference in width of two adjacent steps, l n — l n -\. 
This difference shows a maximum at time t = n. However, the height of the peak decreases 
with increasing n, which means that this effect decreases quantitatively with n. 



V. NUMERICAL RESULTS 

The nonlinear character of the evolution of the crystal shape, as shown by eq.(^) makes 
it impossible to find a general analytic solution of the shape as a function of time. In the 
previous sections we have found such analytic expressions in different asymptotic limits. 
However, in order to study the evolution of the profile at any time, it is necessary to carry 
out a numerical study of the system (p. The results will enable us to check the different 
approximations introduced so far, both in the continuum and in the discrete cases. 

We have solved the differential set of equations ([!]) corresponding to the discrete model 
in the symmetric case, and initial conditions given by eq. (|26|) . We have used a fourth order 



Runge-Kutta algorithm [fTH| . We have considered a set of 200 steps. The neighbouring 
left-side step to our first step is considered to be at infinity. The rightest one is assumed 
to move at constant speed, which is our boundary condition and means that the results 
are obtained in a reference system which moves with the initial slope of the steps. This 
boundary condition means that the terraces at the right of the last one have not changed 
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significantly their width. This will be true only for a certain time scale. As soon as the 
width of the rightest terrace begins to evolve, our solution will be wrong. From the analytic 
results of previous sections, it may be argued that at t 7^ all terraces change their width. 
However, this change is extremely small and, as explained in the previous subsection, only 
the first ones change significantly, so it takes some time until the terraces at the bottom 
change their width appreciably. This fact gives us an easy way to check the validity of the 
solution. It turns out that all relevant features, at least for the first terraces, take place 
during the time scale in which the numerical solutions subject to our boundary condition 
are valid. However, the numerical study at very short times is obscured by the discrete time 
step. 

The differential equations are solved for different values of the initial width, 1(0), as well 
as different time steps, checking the stability of the solution against computational artifacts. 

Figure |3| shows the profile at different times. Due to the fast increase of the width of 
the first steps with respect to the other steps, it is necessary to choose two different length 
scales in the vertical and horizontal axes. In the vertical axis we take the step height as 
the unit length while in the horizontal axis k" 1 is considered the unit length. This is the 
reason why the initial profile appears as an almost vertical straight line. Moreover, 2Ak 2 is 
the unit time throughout the numerical calculations. The curves are plotted every 9 time 
steps, which corresponds to a plot every At = 0.01. 

We have also studied the behaviour of the surface at intermediate times. In this regime 
the numerical solution agrees with the expressions fll7|) in the symmetric model, and ( p5]) in 
the total inverse Schwoebel model, supporting the symmetry argument we have employed 
in order to derive such equations. In Figure f| we show the width of the first three terraces 
as a function of time, and both the logarithmic dependence and the numerical factors are 
recovered exactly after a short transient in both models. The difference in the slope is easily 
appreciated. This good matching with the analytic predictions is again in agreement with 
the idea that the solution ( PD| ) is valid for short times, but that the saturation regime is 
achieved in few time units (in the above mentioned units). 
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We can compare the shape of the discrete model with the predictions of the continuum 
approximation. In Fig. |3] we have checked the validity of Frank's construction for the discrete 
model. We find that all lines joining points of the crystal surface where the slope has a given 
value at different time meet at a point O, which turns out to be the angular point at t — 0. 
Thus all shapes are homothetic of a particular one, and the homothety center is O. 

Finally, in Figure f| the self-similarity of the profile is checked. Neglecting a short tran- 
sient one observes a perfect scaling of the whole profile. The continuum model predicts three 
scaling regions, corresponding to the two initial straight lines and the part of the profile that 
due to the dynamics has deviated significantly from the other two. Therefore, two different 
crossovers are expected. However, they are not readily observed numerically because of the 
finite time step which produces a spurious propagation front which blurres such a crossover. 



In this section we will consider the evaporation of a solid limited by a periodic array of 
grooves (Fig. || a). The sides of the grooves will be assumed to be planar at the initial time 
t = 0. It is sufficient to study the evolution of a half period (Fig. |6] b), which contains a 
time-dependent number n max {t) of steps. The origin x = will be chosen at the highest 
point, supposed to be the left hand side of the half-period. The steps will be labelled 1, 
2,...,n max from the left to the right. 

Only the symmetric case will be considered, so that the evolution of an isolated step is 
governed by equations (|1|) and (|5|). The two lowest steps of each period deserve a special 
attention because they move in opposite directions, so that a facet appears at the lowest 
parts of the profile (Fig. ^| c), Therefore, the position n max of the last step depends on the 
time t. 

It may be of interest to rewrite the equations of motion for the uppest and lowest steps 

as: 



VI. GROOVES 




(34) 
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-An 



tanh ( «(£ + Xnm J\ + tanh ( %n —~ l 



(35) 



where L is the initial width of the groove. Two new features appear with respect to the 
evolution of the profile studied in the previous sections. On one hand, the highest step will 
disappear after some time, and at that time it is necessary to update the x n 's and to replace 
n ma x by n max — 1. A similar updating will take place every time x\ vanishes, until the last 
step disappears and the surface becomes perfectly smooth. On the other hand, since a facet 
appears at the bottom, the first term at the right hand side of equation fl35|) becomes large. 
Therefore, the lowest step moves faster than the other steps and, after some time, it reaches 
the step just above it. Evaporation generates step bunches! Step bunching is frequently 
observed after growing or annealing crystals, and various explanations have been proposed. 
The present one, though very simple, seems to be new. 

When the lowest step has reached the next one, equation ([35D should be modified, because 
the step positions should satisfy the conditions x n > x n -\- The new equation may be found 
if one assumes, following BCF, that steps (forming bunches or not) are in equilibrium with 
the bulk. This results from the assumption that the motion of atoms along steps requires 
lower activation energies than atom detachment from steps. It follows that the chemical 
potential on the lowest terrace (the broad one) near a step should be the bulk chemical 
potential, and that the adatom density on the lowest terrace near a step should be the 
equilibrium adatom density, just as in the case of a single step. On the other hand, the 
adatom density p[x) satisfies the diffusion equation (A.l) (see appendix) independently of 
the number p of steps in the bunch which limits the lowest terrace. Since the boundary 
condition is also independent of p, p(x) is independent of p. Therefore the current at the 
right hand side of the lowest step, which is the gradient of p, is also independent of p, and 
therefore it is tanh(«;(L/2 — x nmax )). On the other hand, the current to the left of the highest 
step in the lower bunch is tax\h(K(x nmax ^ p — x nmax _ p+ i)/2) independently of p. There is a 
current of atoms inside the bunch, the effect of which is to maintain the chemical potential 
uniform within the bunch. It is seen that the current from the lowest step is larger than 
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from the highest one. This ensures the stability of the bunch, i.e., its highest step does not 
move faster than the lowest one. Therefore all velocities within the bunch are the same: 

Ak 
V 



tanb ( «(- - x n )) + tanh ( fi !^=^l!=±l 



(36) 



2 

If bunches of steps appear at other places than the bottom of the profile, equation ([[]) 
should be modified and replaced by an equation similar to (pop. However, no large bunches 
have been observed in the numerical solution of the equations, so that the resulting profile 
is essentially that represented on Fig. || d. If a bunch appears, it is necessary to check its 
stability. As said above, a bunch is stable if the current from its lowest step is larger than 
the current from its highest step. Thus, the condition for a bunch of p steps to be stable if 
its highest step is step m is 

tanh (^ x m-^m-i ^ < _±_ tanh ^ x m+p+1 -x m+p ^ (37) 



We have performed a numerical study of the evolution for the profile shown in Fig. [L0 
a, with 1(0) = 0.1k. Initially, the highest and lowest terraces are twice as broad as the 
others. Fig. depicts the profile at different times during evaporation. The top of the 
profile remains linear while, at the bottom, a bunch appears as expected. At the beginning, 
this bunch grows, but, in our calculations, it reaches a maximal size after some time. This 
size is reached when the width of the lowest terrace becomes of order 1/k so that further 
increase of that width does not produce a further increase of the current from the lowest 
step. Thus, the velocity of the bunch is of order 1/p according to (p6|). This velocity should 
be about the same as that of the next step, which is of order 1. Therefore, the number p 
of steps in the bunch is given by p = l/(/(0)/t). After the bunch has reached its maximal 
size, a pairing of steps is observed. We explain the formation of a step pair just above the 
bottom bunch as follows. Just after reaching saturation, the bunch can have (because of the 
high p value in equation fl36|)) a velocity smaller than the next step. Therefore this step will 
move at a higher velocity than the others and will form a pair of steps with the next one. 
The formation of the other step pairs are presumably of similar origin. Step pairing is not 
always observed and depends on the initial conditions. 
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It is easily deduced from the above arguments that, if the initial width is larger than 1/k 
no bunches are observed. 

It is interesting to compare our results (Fig. |6] and |7|) with what might be expected 
from the continuum approximation. For the sake of simplicity and because it corresponds 
to usual experiments, only the case I « 1/k will be addressed. 

Then, equation ([H]) reduces to 

z = An 2 (38) 

so that evaporation just translates the profile downwards. However, the lowest terrace does 
not evaporate, so that the resulting profile is essentially that of Fig. ^|c. A similar conclusion 
would be obtained from Frank's construction. The continuum approximation fails to predict 
the bunches observed in Figures || d and [7|. As in the previous sections, the failure of the 
continuum approximation is related to the existence of angular points, namely those at the 
bottom of Figure |6| c. However, in the profile studied in the previous sections the initial 
angular point became smooth, and therefore Frank's construction applied at latter times, 
while in the sawtooth profile of Fig. ^ the angular points remain and even give rise to facets 
which are unexpected in the continuum approximation. 

We have also studied the evolution of initially sinusoidal profiles or other differentiable 
profiles (Fig. |] a and |8| b). The main new feature is that an angular point appears at the 
top if the upper terrace is broader than 1/k (Fig. [8] a). The reason is the saturation effect 
associated with the hyperbolic tangent in fl5|): the top of the profile does not evaporate while 



the remainder of the surface follows equation (|3~8). This angular point will nevertheless not 



be observed in most of real materials because large surface adlacunes (not taken into account 
in the present model) will be nucleated. 

We have also made preliminary studies of the evolution of periodic profiles with Schwoebel 
effect. The situation is rather complicated because of instabilities, reminiscent of those found 
by Kandel and Weeks |T!| in a model combining an extreme Schwoebel effect, growth and 



impurities. 
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VII. CONCLUSION 



We have studied the evaporation of a crystal within the Burton-Cabrera-Frank model 
in two cases: i) when the surface is limited by two planes, one of them having a high- 
symmetry orientation (Fig. |l] a); ii) for a "grooved" surface (Fig. ^j. We have paid attention 
to the Schwoebel effect (asymmetry of the sticking coefficient) and to the validity of Frank's 
theorem, based on the continuum approximation. However, we have not investigated cases 
where the instabilities described by Schwoebel || and by Kandel and Weeks [[HJ occur. 

In the case of a corner (Fig. [T] a), Frank's construction can only predict the evolution of 
the crystal above a particular plane II, and below another plane II', but not between the 
two planes. Our numerical solution shows a rounding of the corner and an evolution of the 
profile towards a self-similar shape. We have also shown that the distance between upper 
ledges diverges logarithmically with time except in the case of a total normal Schwoebel 
effect. This divergence wouldn't occur if the continuum approximation were exact. An 
approximate form of the equation of motion has also been investigated, which allows an 
exact solution. 

In the case of a grooved surface, the bottom of the grooves is found to flatten, but their 
edges become steeper due to step bunching. The mechanism responsible for this effect is a 
very simple one, however it cannot be deduced from the continuum approximation. 

The present work can help to understand the formation of defects during annealing of 
crystal surfaces. However, direct comparison with experiment is not possible because, in 
the present work, steps are assumed to be straight, and dislocations are neglected. This is 
only correct for small lengthscales, usually around 0.01 cm or less. We have also ignored the 
effect of vacancies, which are certainly important near the melting point of elements [[J , and 
can even suppress the saturation effect appearing in formula (|5|) when kI > 1. Note that 
this saturation effect has been observed experimentally in some materials, as for example 



did Keller in NaCl |13 |. Neglecting vacancies is only correct if terrace sizes are small. A 
quantitative discussion has been given in Ref. |14| . 
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It is interesting to compare our results with the ones reported by Stoyanov || for a 
stepped surface in which an external force acts on the adatoms. In that case, bunch formation 
is predicted when ^(l) = <p r {l) ~ is positive for a terrace wider than its neighbouring 
terraces, of width I. The sign of depends on the direction of the external force. In our 
case, for an initial surface given by Fig. [I] a, \I/(7) is positive in the case of inverse Schwoebel 
effect. Bunches have not been observed because they will develop above the terrace which 
is larger than its neighbouring ones. In our case, these would correspond to terraces on the 
top of our first terrace and have not been considered in the analysis. What regards bunch 
formation in grooves, it is observed even in the symmetric case, although <p r (l) = <fi(l), 
because \I/(Z) is always positive for the terrace at the bottom. 
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APPENDIX A: VELOCITY OF THE STEPS WITH SCHWOEBEL EFFECT 

In this appendix we derive the expression for the step velocity as sketched in section 1. 
If the diffusion process of adatoms on the surface is much faster than the motion of steps, as 
supposed by BCF 0, the adatoms reach a stationary state during the motion of the steps 
and the velocity of the step is related to the flux of such an adatom density. Then, we should 
determine the density profile on the terraces in order to calculate the step velocity. 

Between two steps, the adatom density p n (x) satisfies the following equation 

at ox 1 r 

x n <C x <C x n j r \ (Al) 

where in the right hand side the first term corresponds to diffusion, the second one to 
evaporation and the third one to deposition. In the present article, F is taken to be zero. 



The boundary conditions at the steps are the following ones, which state the equality of 
two expressions of the current density j(x): 

- j(x) = = D"p - D" Pn (x) (x = x n - e) (A2) 

ox 

-j(x) = = D'p - D'p n {x) (x = x„_! + e) (A3) 

ox 

(A4) 

where e denotes a small quantity. -D'po and D"po are the current of adatoms detaching from 
both sides of the step. These values are imposed by detailed balance. 

Diffusion is usually much faster than step motion, so that the left hand side of (|A1|) may 



be replaced by zero. It is then straightforward to calculate p„(x), and to deduce the current 
density on both sides of each step using eqs. (|A2|) and (|A3|) , which are precisely (pi and (p r . 
Equations @ and (|3|) are therefore easily deduced after some algebra. 

When I is large, the second order expansion of (fj) and (HD in powers of exp(— nl) is given 



by eqs. (|12l) and where 



1 "r" D'D" D' D" 

A ' = 1 + w (A6) 



D' V D' J 1 + K D + «° + 

\ ^ D'D" ' D' D" 

and A" and C" are obtained by interchanging D' and in ( |A6| ) and ( |A7| ) respectively. 

APPENDIX B: SHORT-TIME SYMMETRIC CASE PROFILE 

In this appendix we derive eqs.(|3"U[), starting from the approximate terrace width evo- 
lution equations fl2Tf). Using vector notation these are rewritten in the more convenient 
form 

\l(t) >= ^AK 2 Y\l(t) > +Ak\1 > (Bl) 
with the vector \l(t) > being 

19 



11(f) >- 



( h(t) \ 
hit) 

ln{t) 
\J>N{t) ) 



(B2) 



and where we have introduced 



Y 



-1 .. 

1 0-100.. 
10-10.. 
10-1.. 
1 0.. 



V 



and 



II >= 










(B3) 



The formal solution of ( pi|) is easily checked to be 



\l(t) >= e nt (\l(0) > +2k- 1 Y- 1 \1 >) - -y -1 |l > 



(B4) 



with Q = \Ak 2 Y. 



In order to find an explicit solution to eg. (p3|) , we should decompose the matrix Y in 
eigenvectors, so that the exponential e Qt can be diagonalized. To this end, we introduce the 
vectors 



< k\ 



J-, o , o , . . . , o 



-ink 



iNk 



(B5) 



with n a natural number, n = 0, 1, . . . , N. We can now construct the eigenvectors of Y, 
which is a finite Toeplitz matrix, as an appropriate combination of \k >. It is readily checked 
that 



lib) = (e lk \k > +e- lk \7T -k>) 



(B6) 
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are such vectors if N is an odd number, and k takes the values k = ^(n+i) ' 3 = ~ •••> -W+l- 
The corresponding eigenvalues are — 2zsin(A;). Then, the exponential, in terms of these 
eigenvectors, has the form 

e \A K HY =J2e~ iAK2tBhl{k) \k)(k\ (B7) 
k 

We should now express the exponential in real space. To this end, we have to evaluate 
the action of the operator on vectors |n >, which has its n-th component equal to 1, and 
the rest of them equal to zero. Using the fact that 

< n \k) = -j= (e ink + (-l)*-^-*"*) (B8) 

one can determine the value of the exponential acting on such vectors 

< n\e^ AK2tY \m >= ^ e -^ 2 *-W < n \ k )(k\m > 

k 

\ " iAn 2 t sin(fc) ( '—ink i / i \n—l Ank\ ( Amk i / -i \m— 1 — imk\ lT>c\\ 

= 7^1^ e V + e ) [ e + e ) l By ) 

k 

As we are interested in the situation where there is a large number of terraces, and then 
the behaviour of the system will not be very sensitive on the specific value of N, which is 
large, one can approximate the sums in eq.( P9|) by integrals 



< n\e^ AKHY \m >= — f^ 2 dke- iAKHshl{k) x 

2-K J-tt/2 

^ e i(m-n)k _|_ ^_^(ri+rn) e i{n-m)k _|_ ^_^(m+l) e -i(m+n)k _|_ (n+1) e i(m+n)k^ (BIO) 

Then, using the equality 

— I"' 2 dke- lAft2tsin{k) {e ipk + {-l) p e- ipk ) = J p {-An 2 t) (Bll) 

2lX J-tt/2 

with J n (x) being the Bessel function of order n, one determines the value of eq.( |B10|) , namely 
< n\ exp(AK 2 tY/2)\m >= J„{AkH) + (-l) m+1 J n+m (AK 2 t) (B12) 
which enables us to write down the expression of the exponential in real space 
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exp(An 2 tY/2) 



J0 + J2 —J1 — J3 J2 + J4 — 

J\ + ^3 Jo — J4 —Ji + J5 J2 — J% 

J2 + Ja J\ ~ J5 Jo + J& ~J\ ~ J7 

J3 + J5 J2 — Je J\ + Ji Jo ~ J$ 



(B13) 



v ^ ^ ^ J 

where all the Bessel functions have an argument equal to An 2 t. 

On the other hand, when applying Y on \l(t) > the formal solution 



is expressed 



-h(t) 
h(t)-h(t) 
k{t)-h{t) 



V 



exp(AK 2 tY/2) 








/ 



V 



+ 2kT 1 |^ > 



/ 



(B14) 



/ 



Now, applying expression (|B13| ) for the matrix to eq. ( B14|) we finally arrive at expressions 
(j30l) for the terrace widths as a function of time and the number of the terrace, once the initial 
terrace width is known. 
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FIGURES 

FIG. 1. a)Profile of a crystal surface made of two semi-infinite half planes. The left-side surface 
corresponds to a high symmetry plane, b) Scheme of the movement of adatoms on terraces, showing 
the probability of each event. Schwoebel effect is taken into account through the coefficients D' 
and D". 

FIG. 2. Frank's construction. Left hand partxross section of the initial crystal (dashed curve) 
and of the actual one (full curve); right hand part: the surface generated by the point M defined 
by = n/v(n) . The point P of the crystal surface where the normal is n moves on the dotted 
straight line, parallel to the normal to T at M. On the left hand of the dotted line, the crystal 
is planar and parallel to its original orientation. Only the useful part of the curve of M has been 
shown. Note the presence of a point at infinity. It corresponds to I — ► oo and its abcissa is 
— 1, corresponding to the formula QM = —(l,l)/ip(l), where tp is given by (5). 

FIG. 3. Snapshots of the evolution of the initial profile shown in Fig. 1 a in the symmetric 
case. Profiles are drawn every 20 times in the units defined in the text. Two different scales have 
been chosen in the explained in the text. The profile in the discrete symmetric case is 

shown to be in agreement with Frank's construction. Two equal slopes at different times are seen 
to be joined by a straight line, and different straight lines converge at 0, the angular point at t=0. 

FIG. 4. Width of the three first terraces as a function of time at long times. A, o and y 
correspond to the symmetric case, while o, • and □ correspond to the total inverse Schwoebel 
effect. The logarithmic behavior, the prefactor and the constant are recovered in both cases. The 
1/2 difference in the slope is easily appreciated. 

FIG. 5. Self-similarity of the profiles of the crystal without Schwoebel effect. The crossovers 
between the different regimes are blurred by numerical artifacts. 
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FIG. 6. a) Initial profile of the surface studied in Section 6. b) Enlargement of the half-period 
inside the rectangle of Fig. 6 a. c) The evaporation shape as would be predicted from the continuum 
approximation, d) The real shape resulting from step bunching, not taking into account the step 
pairing occasionally observed in simulations. 

FIG. 7. Snapshots of the evolution of the initial profile shown in Fig. la. Profiles are drawn 
every ten times in the units defined in the text. Two different scales have been chosen in the axes, 
as explained in the text. 

FIG. 8. Snapshots of the evolution of an initial parabolic profile, a) Evolution when initially 
all terraces are smaller than 1/k. b) Evolution when initially the top terraces are larger than 1/k 
and the bottom ones smaller than this value. 
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